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Abstract 

We present two methods for computing dimensionaily-reguiated NRQCD heavy-quarkonium 

m t 

! matrix elements that are related to the second derivative of the heavy-quarkonium wave function 

at the origin. The first method makes use of a hard-cutoff regulator as an intermediate step 



m 

O 
O 



and requires knowledge only of the heavy-quarkonium wave function. It involves a significant 
cancellation that is an obstacle to achieving high numerical accuracy. The second method is more 



direct and yields a result that is identical to the Gremm-Kapustin relation, but it is limited to use 
in potential models. It can be generalized to the computation of matrix elements of higher order in 

& ■ 

the heavy-quark velocity and can be used to resum the contributions to decay and production rates 
that are associated with those matrix elements. We apply these methods to the Cornell potential 

X 

model and compute a matrix element for the J/ip state that appears in the leading relativistic 
correction to the production and decay of that state through the color-singlet quark-antiquark 
channel. 
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I. INTRODUCTION 



In the effective field theory nonrelativistic quantum chromodynamics (NRQCD), the lead- 
ing relativistic corrections to S-wave heavy-quarkonium decay and production processes in 
the color-singlet quark-antiquark channel are proportional to matrix elements that are re- 
lated to the second derivative of the quarkonium wave function at the origin. These matrix 
elements are inherently nonperturbative in nature. Their importance in phenomenological 
calculations has led to a number of attempts to determine their values. 

Even before the introduction of the NRQCD approach for quarkonium decay and produc- 
tion [1] , these matrix elements appeared in phenomenological studies of charmonium decays 
[2]. Owing to uncertainties that arise from the uncertainty in the charm-quark mass m c 
and from uncalculated terms of higher order in the quantum chromodynamic (QCD) strong 
coupling a s , such phenomenological determinations have not led to accurate values for the 
matrix elements. There also have been attempts to determine the matrix elements in lattice 
calculations [3]. In this case, large uncertainties arise because there is a substantial can- 
cellation that occurs when one converts from lattice to continuum dimensionally-regulated 
matrix elements. In principle, one can determine these matrix elements approximately by 
making use of the Gremm-Kapustin relation [4], which expresses the matrix elements in 
terms of the quarkonium and heavy-quark masses. (See Ref. [5] for an example of this ap- 
proach.) Unfortunately, this method is plagued by large uncertainties in m c . In both the 
lattice and Gremm-Kapustin approaches, the uncertainties are so large that even the signs 
of the matrix elements are in some doubt. 

A further difficulty that complicates the calculation of the matrix elements that are 
related to the second derivative of the wave function at the origin is that they contain a linear 
ultraviolet (UV) divergence, and, hence, must be regulated. Dimensional regularization of 
these matrix elements is particularly useful because it is consistent with existing calculations 
of quarkonium decay and production rates at relative order a s and a 2 s . 

In this paper, we present two methods for calculating the second derivative of the wave 
function at the origin. In the first method, we initially use a hard-cutoff regulator, with 
cutoff A, to define the relevant matrix element. Then we compute the difference between 
the hard-cutoff regularization and dimensional regularization in perturbation theory. We 
subtract this difference from the hard-cutoff matrix element. There remains a dependence 
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of the matrix element on A that falls as 1/A in the limit A — > oo. That dependence can be 
removed by calculating at a number of values of A and extrapolating to A = oo. However, 
in this method, it is difficult to achieve high numerical accuracy because large, cancelling 
A-dependent contributions appear. The second method that we present bypasses the hard- 
cutoff step, but it is applicable only to potential models for the wave function. It yields a 
result that is identical to the Gremm-Kapustin relation [4] for the potential model. That 
result can be used to resum certain contributions of higher order in v to amplitudes that are 
computed in NRQCD. 

Having established a formal procedure for computing the relevant matrix elements, we 
carry out an explicit computation for the J/ip (or i] c ) states in the Cornell potential model. 
We do not distinguish between the J/x/j and rj c wave functions, which differ only in corrections 
of relative order v 2 , where v is the velocity of the heavy quark or antiquark in the quarkonium. 
(v 2 ~ 0.3 in charmonium and v 2 ~ 0.1 in bottomonium.) In principle, if we know the static 
heavy- quark- antiquark {QQ) potential exactly, then we can calculate the quarkonium wave 
function of the leading QQ quarkonium Fock state up to corrections of relative order v 2 . 
Existing lattice data for the static QQ potential yield values for the string tension. We 
examine values for the parameters in the Cornell potential that bracket the lattice values 
for the string tension. 

In our numerical calculations, the results from our two approaches agree well and give 
a value for the second derivative of the wave function at the origin that is in agreement 
with expectations from the f-scaling rules of NRQCD [1]. The largest uncertainties in our 
calculation are of relative order v 2 . Therefore, our determination of the second derivative 
of the wave function at the origin is the most accurate to date and should be useful for 
phenomenological studies of quarkonium production and decay. 

The remainder of this paper is organized as follows: In Sec. II, we give a brief description 
of the Cornell potential model. In Sec. Ill, we discuss the NRQCD matrix elements that 
are relevant to this work. Section IV contains a description of the hard-cutoff regulator. In 
Sec. V, we explain how we compute the difference between hard-cutoff regularization and di- 
mensional regularization. The direct method of calculation, which bypasses the hard-cutoff 
step, is discussed in Sec. VII, along with its application to resummation of contributions of 
higher order in v. In Sec. VI, we decompose the difference between a hard-cutoff matrix ele- 
ment and a dimensionally-regulated matrix element into sums of short-distance coefficients 
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times dimensionally-regulated matrix elements. Such a decomposition is used in existing cal- 
culations of the difference between a lattice matrix element and a dimensionally-regulated 
matrix element. Sec. VIII contains our numerical results and a discussion of them. We 
summarize our results in Sec. IX. In Appendix A, we illustrate our methods for the case of 
a pure Coulomb potential, and, in Appendix B, we list some integrals that are useful in our 
analyses. 

II. POTENTIAL MODEL 

We compute matrix elements for the J/ip or r\ c states using a potential model. In this 
model, we neglect the effects of the heavy-quark spin, which are suppressed as v 2 . There- 
fore, we do not distinguish between the J/ip and r] c wave functions or matrix elements. 
We note that, if we knew the heavy-quark potential exactly, then we could calculate the 
heavy-quarkonium wave function in a potential model up to corrections of relative order v 2 
(Ref. [6]). We make use of the Cornell potential model of Ref. [7]. For appropriate choices of 
parameters, the Cornell potential provides a reasonably good fit to heavy-quark potentials 
that are measured in lattice calculations. 1 

Now we summarize the methods that we use to constrain the parameters of the Cornell 
potential and to solve the Schrodinger equation. We refer the reader to Ref. [9] for further 
details. 

The Cornell potential [7] is given by 



where the parameters k and a determine the strength of Coulomb and linear potentials, re- 
spectively. For a color-singlet QQ pair, the Coulomb-strength parameter k can be expressed 
in terms of an effective strong coupling a s as 



where Cp = 4/3. The parameter a is related to the string tension a as 

1 

1 For a recent review that discusses heavy-quark potentials from lattice measurements, see Ref. [8]. 



V(r) = -- + - 2 , 
r a 1 



(1) 



k = a s Cp, 



(2) 
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Following Ref. [7], we replace k and the binding energy ee by dimensionless parameters 
A and (: 

2 

k = (ma)" A, (4a) 

4 

e B = m(ma)"3(. (4b) 

For a given value of A, we fix the heavy-quark mass m and the parameter a in the Cornell 
potential as follows. First, we require that the energies that result from the solutions to 
the Schrodinger equation match the measured values of the difference of the J/ip and ^>(2S) 
masses. We use Mj/^ = 3.096916 GeV and M^s) — 3.686093 GeV. Second, we require that 
the wave function at the origin matches a value that is derived from the measured value of 
the leptonic width of the J/ip and the perturbative formula 

r[J/ ^ rt -] = ^ W 0)p(l-^) 2 . (5) 

Here, tp(0) is the wave function at the origin and e c = 2/3 is the fractional electric charge of 
the charm quark. In Ref. [5], ip(0) is estimated by using the the formula (5) at both leading 
order in a s (LO) and next-to- leading order in a s (NLO). The results are 

f 0.18619 GeV 3/2 (LO) 
^(0) = I K ' (6) 

y 0.23629 GeV 3/2 (NLO), 

where, for convenience, we have taken ip(0) to be positive and real. In order to estimate the 
effects of the uncertainty in ip(0), we carry out our calculations for both the LO and NLO 
values of ip(0). For a more detailed discussion of the determination of m and a, see Ref. [9]. 

Values of m, a, and the scaled energies ( of the IS and 2S states for various values of A 
are shown in Tables I and II, along with values for a s from Eq. (2), a from Eq. (3), and 7c 
from Eq. (12). Table I contains values of the potential-model parameters that correspond 
to the LO value of ip(0), while Table II contains those that correspond to the NLO value of 

Lattice measurements of the heavy-quark potential yield values for effective coupling 
a s of 0.22 in the quenched case and approximately 0.26 in the unquenched case [8]. A 
lattice measurement of the string tension K = a (Ref. [10]) gives Kal = 0.0114(2) at 
a lattice coupling (5 = 6.5, where is the lattice spacing. Lattice calculations of the 
hadron spectrum at (3 = 6.5 yield values for l/a L of 3.962(127) GeV (Refs. [11, 12]) and 



TABLE I: Potential-model parameters and derived quantities as a function of the strength A of the 
Coulomb potential. The definitions of the parameters and derived quantities are given in the text. 
The parameters are computed using the inputs Mj^ = 3.096916 GeV, M^s) = 3.686093 GeV, 
and the LO value V(0) = 0.18619 GeV 3 / 2 , as is described in the text. 



A 





0.2 


0.4 


0.6 


0.7 


0.8 


1.0 


1.2 


1.4 


ClO 


2.33811 


2.16732 


1.98850 


1.80107 


1.70394 


1.60441 


1.39788 


1.18084 


0.95264 


C20 


4.08790 


3.97017 


3.85003 


3.72747 


3.66528 


3.60249 


3.47510 


3.34529 


3.21307 


m (GeV) 


1.70670 


1.51548 


1.35120 


1.21003 


1.14710 


1.08877 


0.98458 


0.89501 


0.81796 


a (GeV" 1 ) 


1.97932 


2.08520 


2.19805 


2.31833 


2.38139 


2.44648 


2.58295 


2.72816 


2 88253 


(GeV) 


0.50522 


0.47957 


0.45495 


0.43134 


0.41992 


0.40875 


0.38715 


0.36655 


0.34692 


a (GeV 2 ) 


0.25525 


0.22999 


0.20698 


0.18606 


0.17633 


0.16708 


0.14989 


0.13436 


0.12035 




0. 


0.06966 


0.14519 


0.22624 


0.26866 


0.31225 


0.40255 


0.49634 


0.59272 


7c (GeV) 


0. 


0.07037 


0.13079 


0.18250 


0.20545 


0.22664 


0.26423 


0.29615 


0.32322 


TABLE II: As in Table I, except that the potential-model parameters are computed using 


the NLO 


value ip(0) - 


= 0.23629 GeV 3 / 2 . 
















A 


0.2 


0.4 


0.6 


0.8 


0.9 


1.0 


1.1 


1.2 


1.4 


C10 


2.16732 


1.98850 


1.80107 


1.60441 


1.50242 


1.39788 


1.29071 


1.18084 


0.95264 


C20 


3.97017 


3.85003 


3.72747 


3.60249 


3.53910 


3.47510 


3.41050 


3.34529 


3.21307 


m (GeV) 


2.08228 


1.85655 


1.66259 


1.49597 


1.42168 


1.35282 


1.28896 


1.22975 


1.12388 


a (GeV" 1 ) 


1.92598 


2.03021 


2.14130 


2.25967 


2.32171 


2.38571 


2.45174 


2.51984 


2.66243 


(GeV) 


0.51922 


0.49256 


0.46700 


0.44254 


0.43072 


0.41916 


0.40787 


0.39685 


0.37560 


a (GeV 2 ) 


0.26959 


0.24262 


0.21809 


0.19584 


0.18552 


0.17570 


0.16636 


0.15749 


0.14107 


a s 


0.05942 


0.12387 


0.19301 


0.26638 


0.30448 


0.34342 


0.38310 


0.42343 


0.50566 


7c (GeV) 


0.08249 


0.15331 


0.21393 


0.26567 


0.28859 


0.30972 


0.32920 


0.34714 


0.37887 



3.811(59) GeV (Refs. [11, 13]). These yield values of the string tension of K = 0.1790±0.0119 
and K = 0.1656 ± 0.0059, respectively. 

Comparing the results of these lattice measurements with the LO parameters in Table I, 
we see that the values of the string tension at A = 0.7 and 0.8 span the range of lattice results 
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for the string tension, while the values of a s at A = 0.6 and 0.7 span the range of lattice 
results for a s . Comparing the results of the lattice measurements with the NLO parameters 
in Table II, we see that the values of the string tension at A = 1.0 and 1.1 span the range of 
lattice results for the string tension. However, the values of a s at A = 0.9, 1.0, and 1.1 are 
all larger than the lattice values. It is not clear whether this discrepancy between the lattice 
and NLO potential-model values for a s arises from the use of an inaccurate value for ip(0), 
from effects due to the running of a s , which are not taken into account in the fits to the 
lattice data, or from the absence of corrections of relative order v 2 in the Cornell potential 
model. However, we note that the NLO values for a s at A = 0.9, 1.0, and 1.1 do not differ 
greatly from the value of the running a s at the scale of the heavy-quark momentum m c v . 

Finally, we mention that we obtain the J/vp wave function by expressing the radial part 
of the Schrodinger equation as a difference equation, which we integrate numerically. See 
Ref. [9] for details. 

III. NRQCD MATRIX ELEMENTS 

In this section we describe the NRQCD matrix elements that are relevant to our calcu- 
lation. 

In the rest frame of an S— wave heavy quarkonium H in a spin-singlet ( 1 5'o) or spin-triplet 
( 3 Si) state, one can express the wave function at the origin of the leading QQ Fock state in 
terms of the following color-singlet NRQCD matrix elements [1]: 



Here ip and x are Pauli spinor fields that annihilate a quark and an antiquark, respectively. 
The bilinear operators involving ^ and x are evaluated at zero space-time position, cr is 
a Pauli matrix, and e is the quarkonium polarization vector. ip(p) is the momentum-space 
wave function for the leading Q{p)Q{—p) Fock state in the quarkonium. The wave function 
is, of course, gauge dependent. Throughout this paper, we work in the Coulomb gauge. 
The normalization factor l/y/2N c accounts for the traces in the SU(2)-spin and SU(3)-color 
spaces. In equating the wave functions for the spin-singlet and spin-triplet cases, we are 
ignoring effects of relative order v 2 . 




(7b) 



(7a) 
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Relativistic corrections to the production and decay rates for a heavy quarkonium involve 
matrix elements that are related to the second derivative of the wave function at the origin: 



^ (2) (0) = J ^P^(p) = 7^°^ (- V2 ) tWSo)), (8a) 

^ (2) (0) = e J -^P^ip) = T^ ^ (" V2 ) ^W^))- (8b) 

Usually, these operator matrix elements are written in terms of the covariant derivative D 
(Ref. [1]), rather than V. However, in the Coulomb gauge, the difference between the D 
and V is suppressed as v (Ref. [1]). 

The quantity ^^(O) is the focus of this paper. It is common in phenomenology to make 
use of a parameter 



V) = A (9) 
{ ' m^(0)' { ' 

where m c is the charm-quark pole mass, which we distinguish from the parameter m that 
appears in our potential model. Note that ^ 2 \0) is different from the expectation value of 

^ {2) (o)^ J ^pV(p)^). (io) 

Let us investigate the ultraviolet behavior of the matrix elements in Eq. (8). At large mo- 
mentum \p\, the interaction of the QQ pair in QCD is dominated by the Coulomb potential. 
In this limit, the bound-state wave function approaches the pure Coulomb wave function: 

Here 7c is a parameter that is related to the binding energy of Coulomb interaction: 



7c = ^a s C F m, (12) 

where a s is the effective strong coupling, Cf = 4/3, and m is the quark mass. Substituting 
Eq. (11) into Eq. (8), we see that the matrix elements in Eq. (8) have a linear ultraviolet 
divergence. Hence, in order to define them, we must impose a regulator. 



IV. HARD-CUTOFF REGULATOR 

Our ultimate goal is to regulate the matrix elements in Eq. (8) using dimensional regular- 
ization. However, as an intermediate step, we impose a hard-cutoff regulator. In principle, 
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the methods that we use could be employed in lattice calculations of the matrix elements, 
in which the lattice supplies the hard-cutoff regulator. However, for the purposes of the 
present work, we make use of a simple, analytic UV regulator in momentum space: 

^=Iwr^"^ < 13 > 

where A is the cutoff. 

We introduce the Fourier transform of the wave function to coordinate space: 

= J ^^P)- ( 14 ) 



For S-wave states, we can write ip(x) = R{r)/y4^r, where r = \x\ and R(r) is the radial 
wave function. Substituting Eq. (14) into Eq. (13), we can carry out the angular integration 
over p to obtain an expression in coordinate space: 



■4 2, (o)= v 



POD 

R(0) - A 2 rR(r)e- Ar dr 
Jo 



(15) 



V. DIFFERENCE BETWEEN HARD-CUTOFF AND DIMENSIONAL REGU- 
LARIZATION 

Now we work out the difference between the hard-cutoff matrix element and the 
dimensionally-regulated matrix element 

A^ 2 )(0)=^i 2) (0)-^ iA (0). (16) 

Because A^^(O) involves a change in the ultraviolet cutoff, it is sensitive only to the high- 
momentum part of the p integration in the momentum-space definition of ^ (2) (0) in Eq. (8). 
Therefore, we can compute Aip^(0) in perturbation theory. Here, we carry out the com- 
putation in lowest (one-loop) order. 2 We note, however, that the perturbative calculation 
of A^ 2 )(0) does not contain all of the corrections that behave as powers of 1/A to the 



2 One could, in principle, compute corrections of higher order in a s . Ultimately, the series of these cor- 
rections diverges, owing to the renormalon ambiguity that appears in dimensionally-regulated matrix 
elements. That ambiguity is canceled by a corresponding ambiguity in the NRQCD short-distance co- 
efficients, provided that one computes to the same loop order in both the matrix elements and the 
short-distance coefficients. See Ref. [f 4] for a discussion of this point. 
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FIG. 1: Feynman diagram corresponding to I^(p), which is the one-loop renormalization of the 
operators in Eq. (8). The solid lines represent the heavy quark and antiquark, and the dotted line 
represents the potential between them. 

(2) (2) 

difference between -0^ (0) and ipj^R aW- Therefore, it is necessary to take the limit A — > oo 
in order to remove the effects of such power corrections. 

We begin by using the Bethe-Salpeter equation (or, equivalently, the Schrodinger equa- 
tion) to expose one explicit loop in the wave function. Then, we have 

^ (2) (0) = / -0f !i2) (P)$(P)i (17) 

where I^ 2 '(p) is the quantity that is represented by the Feynman diagram in Fig. 1. In 
Fig. 1, the dotted line represents the potential. Since the loop integral is dominated by large 
momenta, we can approximate the potential by the Coulomb part. 

We emphasize that we can always make this approximation in QCD, irrespective of the 
use of a potential model. Asymptotic freedom allows one to evaluate the high-momentum 
loop in Fig. 1 in perturbation theory. Then, in the Coulomb gauge, the Coulomb-gluon 
interactions with the heavy quark and antiquark give the result of leading order in v, while 
the transverse-gluon interactions are suppressed as v 2 (Ref. [1]). In the language of NRQCD, 
I^ 2 \p) is the perturbative one-loop renormalization of the operators in Eq. (8) in leading 
order in ot s and v. 

Note that in (p) we retain the binding energy of the heavy quarkonium 

e B = -ll/m. (18) 
10 



We assign momenta such that the Q and Q each carry half of the binding energy in the rest 
frame of the heavy quarkonium. 
is given by 

A^(0) = J ^AJ< 2 >(p)£(p), (19) 

where 

AP\p) = ^\p)-lSl(p). (20) 

Writing out the expression for the diagram in Fig. 1 for the cases of a hard cutoff and 
dimensional regularization, we obtain 

d A k 



A4 2) (p) = A-Ka s C F 



I 



(27T)" 



: (k + P y 



k 2 



An 



2k 

2m 



X 



-k 



7l 



(fc+p) 2 



+ ie 



(fc+p) 2 

2m 

A 2 



+ ie 



(k + p) 2 + A 2 



- 1 



DR 



(21) 



2m 2m 

The first term in the numerator brackets in Eq. (21) corresponds to the hard cutoff, while the 
second term "Idr" corresponds to dimensional regularization. The meaning of the "Idr" 
is that it is unity unless it multiplies a scaleless integral (which vanishes in dimensional 
regularization), in which case it is zero. The subscript "B" in A7g (p) indicates that we 
have retained the dependence on the binding energy — 7 B /2m. Since the integral in Eq. (21) 
is dominated by large momenta, we expect the final result to be insensitive to 7b- However, in 
discussions that occur later in this paper, it is illuminating to retain the full 7b dependence. 

The integral over the loop energy k in Eq. (21) can be carried out by using the residue 
theorem. The result is 



A/ B 2) (p) = Aita s C F m 



d 3 k 



(fc+p) 



A 5 



(fc + p) 2 + A 2 



l-DR 



(22) 



(27r)3fc 2 [(fc+p) 2 + 7 2 ] 

(Note that we also could have reached this expression directly from the momentum-space 
Schrodinger equation.) In the second term in brackets in Eq. (22), we discard a scaleless 
integral of the form 

(23) 



DR 



(2tt) 3 fc 2 ' 

In the first term in brackets in Eq. (22), we rewrite the numerator (fc + p) 2 as [(fc + p) 2 + 
7b] — 7b an d partial-fraction the —7! term. The result of these manipulations is 



A/ B 2) (P) = 



ArcasCpTn C d 3 k 1 
A 2 - 7 2 J (Xjsfc 1 



A 4 



7b 



(fc+p) 2 + A 2 (fc + p) 2 + 7l 



(24) 
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Evaluation of the integral in Eq. (24) is straightforward. The result is 



ai b (p) ~ 



V \ .4 . ( \P\ 



A arctan — — — 7 R arctan . 

A J IB V 7b 



(25) 



|P|(A 2 -7I) 

As expected, I^\p) is insensitive to 7b- Neglecting terms of higher order in 7I/A 2 , we 
obtain 

Aigl(p) = limA/f»( P ) = -^A 2 arctan ( , (26) 

where the subscript "NB" indicates that we have neglected the binding energy by dropping 
contributions of higher order in 7j|/A 2 . (Note that arctan(|p|/7 B ) is well behaved as 7 B — > 
and bounded over the entire range of its argument.) In our numerical analyses, we make 
use of AJnb(p), rather than AI^\p), and consistently neglect the binding energy in short- 
distance (high-momentum) quantities. Of course, the binding-energy dependence in the 
wave function is retained fully. 

In our numerical analyses, we solve the Schrodinger equation in coordinate space rather 
than in momentum space. The Fourier transformation of the coordinate-space wave function 
to momentum space involves an oscillating integrand and, hence, is difficult to evaluate 
numerically. Therefore, it is convenient to evaluate A^^(O) in coordinate space: 

A^ (2) (0) = J d 3 xAlW(x)4>(x), (27) 

where 

A/ (2) (*) = J ^ e *-A/< 2 >(p) (28) 

is the Fourier transform of AI^(p). It is a simple matter to evaluate A/( 2 )(a:) analytically. 
The results are 

r(2)^ _ f Jyx A m ( x _ a s C F m / A 4-Ar 



AI * {X) = J AI * (P) = 4^(V- 7 |) ( e - -Ik- 1 ") > (2*0 

Substituting A/( 2) (a;) into Eq. (27), we obtain 

M 2> (0) = [ (AV* - ifie-"') ^ *', (30a) 



OG 



A^l(0) = a s C F mA 2 I e ~ Ar ^ldr. (30b) 

V47T 
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Then, we obtain the dimensionally-regulated matrix element by computing 

^1a(0)=^1 2) (0)-A^(0) (31) 

and taking the limit A — > oo in order to remove uncompensated power corrections that 
vanish as 1/A: 

^gi(O) = lim ^ >A (0). (32) 

VI. DECOMPOSITION OF A^( 2 >(0) 

It is illuminating to decompose A^ 2 \Q) into products of NRQCD matrix elements times 
short-distance coefficients. In existing lattice calculations [3], the difference between matrix 
elements in hard-cutoff regularization and those in dimensional regularization is expressed in 
terms of such a decomposition. We decompose Aip( 2 \0) by applying the method of regions 
[15] to the integration over the wave-function momentum p in Eq. (19). We take for AI^\p) 
the no-binding expression AI^(p), which is given in Eq. (26). 

First, we extract the leading term, which is obtained from the \p\° part of AI^l(p), 
namely, agCpTnA. It yields a contribution to Aip^(0) that is 

A^ 2) (0) = a s C F mA J ^^(p) = a s C F mA^(0) = 2 7c A^(0). (33) 

We identify this as a one-loop short-distance coefficient times the matrix element of leading 
order in v, namely, ip(0). A^^(O) is linearly divergent in the limit A — > oo and cancels the 
linear divergence in the hard-cutoff matrix element resulting in a matrix element 

■0dr(O)) which, according to the NRQCD f-scaling rules, is of order (mt>) 2 -?/>(0). The v 
scaling of this quantity is verified by the result of an analytic calculation of V'dr(O) f° r the 
case of a pure Coulomb potential in Appendix A. 

Next we examine the remainder of AI^(p), namely, a s CjrmA[(A/|p|) arctan(|p|/A) — 1]. 
We decompose the loop integration over p into regions of small |p|, in which \p\ <C A, and 
large |p|, in which |p| ~ A. 

In the small- |p| region, we expand in powers of \p\. The resulting contributions, which 
we call A^4^(0), have the form 

A4 2 \0) = a s C F m I 7^(p)f>n^£, (34) 

JDR W n=3 A 
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where c n is the coefficient of (|p|/A) n in the power-series expansion of arctan(|p|/A), and 
the subscript DR on the integral indicates that any UV divergence is to be regulated dimen- 
sionally. Aip^\o) corresponds to one-loop short- distance coefficients times dimensionally- 
regulated NRQCD matrix elements of higher order in v. The contributions to Aijj^\o) are 
of order a s m{mv) n ^ 1 / A n ^ 2 and, hence, are suppressed relative to ^dr(O) by at least one 
power of m/A as A — > oo. 

In the large- |p| region, we approximate ip(p) by its asymptotic form at large \p\, which 
amounts to neglecting corrections of higher order in a s and mv/A. The asymptotic form 
of ip{p) is obtained by using the Bethe-Salpeter equation to expose an explicit loop in the 
wave function, as in the derivation of Eq. (22). The result is 

7/ \ f d 3 k ~ 8tt7 C 8tt7 C ^(0) ~ 

*&») ~ J (2^^J^k^ ~ = *»<!»>■ (35) 

Then, we have for the contribution in the large- \p\ region 

A^3 2) (0) = a s C F mA J ^-^^ asy (p) ^ arctan (^^j ~ 

= -(l/2)(a s C F m) 2 ^(0) = -2 7 ^(0). (36) 

A^3 2 ^(0) has the interpretation of a two-loop short-distance coefficient times the lowest-order 
NRQCD matrix element. 

(2) (2) (2) 

Although the contribution Aip^ (0) is suppressed as a s /A relative to A^f'(0), A^ 3 ; (0) 

(2) (2) 

can be important numerically because, in ^ DR (0), Aip{ (0) is canceled by contributions 
from i/j^\o), leaving a small remainder. As we have remarked, V'dr(O) itself is nominally 
of order (mt>) 2 -?/>(0). Therefore, A^^O) is suppressed only as a 2 /v 2 relative to ^^(O). 
According to the w-scaling rules of NRQCD [1], a s is of order v 2 , and, hence, there is a slight 

(2) 

suppression. It turns out that, for the Cornell potential, the numerical value of Aip 3 (0) 

(2) 

is about —40% of Vt>r(0)- This suggests that two-loop corrections to the short-distance 
coefficients could be important numerically in converting results of lattice calculations of 
■?/>( 2 )(0) to continuum regularization. In the case of a pure Coulomb potential, a s is of 
order v, and, so, there is no suppression of A^^O) relative to V'dr(O) at all. We show in 
Appendix A, V'dr(O) satisfies the Gremm-Kapustin relation in the case of a pure Coulomb 

(2) (2) 

potential, while V4)r(0) ~~ ^3 (0) does not. 



14 



VII. DIRECT METHOD OF CALCULATION OF V>dr(°) 

In the method of calculation that we have outlined, it is necessary to take the limit of 
the quantity ^dr,a(0) as A S oes to infinity. V'dr.aW consists of a difference of A-dependent 
terms that grow approximately linearly with A. If one computes at a large enough value of A 
to be near the asymptotic value of ip^R a(0)> then there is a substantial cancellation between 

(2) (2) 

these terms. For example, for the NLO parameters at A = 1.0, the ratio ip A (0)/V4)ra(0) 
takes on the values 18, 34, and 50 when A/m equals 10, 20, and 30, respectively. In numerical 
calculations, the cancellation in ^dra(O) * s a significant obstacle to achieving good accuracy. 
Therefore, it is desirable to have a direct method of computation of V'dr(O) that does not 
pass through an intermediate hard-cutoff step. In this section, we present such a method. 
We also show how it leads to a simple formula for the resummation of a class of corrections 
of higher order in v. 

A. Method 

We begin again by exposing one loop in the matrix element, as in Fig. 1. Now, however, 
since we are computing the matrix element itself, not a difference of matrix elements for 
different regulators, the loop in Fig. 1 is not dominated by large momenta. Therefore, we 
must retain the complete Cornell potential in the corresponding expression. Repeating the 
steps that lead to Eq. (22), but for the complete Cornell potential, and using Eq. (17), we 
obtain for the dimensionally-regulated matrix element 



It turns out that, in the case of the Cornell potential, the binding energy is positive. There- 
fore, it is convenient to define a quantity 7b as 



From the ie prescription in Eq. (37), we see that we can analytically continue from 7b 
positive to 7b positive imaginary. Hence, we make the identification 7b = vjb, where 7b is 
real and positive. Discarding the contribution of the scaleless integral 




(37) 



= -i = S>o. 

m m 



(38) 




(39) 



15 



in Eq. (37), we have 

m7B 7 (27T)3 7 (27T)3 fc2_^2 +2e • l 4U J 

Now we note that ip(p) satisfies the momentum-space Schrodinger equation 

(fc 2 - 7l)^(fc) = -J ^~ 3 mV(k - pty(p). (41) 
Therefore, Eq. (40) can be written as 

^lI(O) = 7l / = 7^(0). (42) 

The Gremm-Kapustin relation [4] follows from the NRQCD equations of motion at leading 
order nontrivial in v. It states that 

(v 2 ) = ^ = - + 0(v% (43) 

Using the definition of 7I in Eq. (38), we see that our result in Eq. (42) is precisely the 
Gremm-Kapustin relation for the potential model with mass m = m c . This is not surprising, 
since a potential model can be derived from NRQCD at leading order in v (Ref. [6]) and, 
hence, implicitly respects the NRQCD equations of motion at leading order in v. 



B. Resummation 

We note that the result in Eq. (42) can easily be generalized to the case of NRQCD 
operator matrix elements involving additional powers of V 2 , such as 

^ (2n) (0) = J -^P^P) = 7^°^ (-V 2 ) n ^(%)>, (44a) 
^ (2n) (0) = e J ^ 3 P 2n kp) = "yf^ ^ (-V 2 ) n Vl#( 3 Si)>, (44b) 

where n is an integer. As we have remarked earlier, such operator matrix elements are 
usually written in terms of the covariant derivative D (Ref. [1]), rather than V. However, 
in the Coulomb gauge, the difference between D and V is suppressed as v (Ref. [1]). 
Repeating the steps that lead to Eq. (37), we obtain 

*w = - m I Swf' {k - (45) 
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In Eq. (45), we write 



JU2n ~2n 

-k 2n - 2 +^k 2n ~ 4 + ■■■ + —% (46) 



k 2 -% + ie m k 2 -% + ie 

and discard scaleless integrals of the form 

/ (0^ (fc - p)fc2n ' lDR = / i^nm+pr^ m 

where n' is an integer. The result is 

miB J (2tt)3 J (2vr)3 fc2-7| + ie W 

where we have used the momentum-space Schrodinger equation [Eq. (41)]. The expression 
in Eq. (48) can be thought of as a generalized Gremm-Kapustin relation that holds in a 
potential model. 

Now, suppose that there are contributions to a decay or production amplitude, computed 
in NRQCD, of the form 

A (2n) = J_#(2n)( ^(2n) (4g) 
71/. 

Here, (l/n\)H^ 2n \0) is an NRQCD short-distance coefficient that is defined by 



(50) 

p=0 



where H^ 2n \p 2 ) is the hard-scattering amplitude for the process, evaluated in the quarko- 
nium rest frame at Q or Q momentum squared p 2 . We can use Eq. (48) to resum the 
contributions A^: 

E A(2n) = E V (2n) ( )*°(0) = E Wb>(0) = # (tIMO). (51) 

* — ' * — ' n! z — ' n! 

n n n 

VIII. NUMERICAL RESULTS AND DISCUSSION 

We now use the hard-cutoff method of Sees. IV and V and the direct method of Sec. VII 
to compute ^dr(O). 

We begin with the hard-cutoff method. We substitute the lS'-state coordinate-space 
radial wave function that we obtain by integrating the Schrodinger equation numerically 
into Eqs. (15) and (30b) and carry out the integrations over r numerically. We then use 
Eq. (31) to compute V'draCO)- The results are shown in Fig. 2. As can be seen, ^\o) grows 
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to the LO potential-model parameters of Table I for A = 0.6, 0.7, and 0.8. The right figure 
corresponds to the NLO potential-model parameters of Table II for A = 0.9, 1.0, and 1.1. In each 

(2) 

figure, the three curves that are nearly linear are 0.1 x ^ ; (0), and the three curves that reach a 

(2) 

plateau at large A/m are j/j^ a (0). 



nearly linearly with A, as is expected from the linear ultraviolet divergence that it contains. 

(2) f2) f2l 

Vt>ra(0) is the result of a substantial cancellation between ip ( A '(0) and A^g(0), at the level 
of about one part in 50 at A/m = 30 for the NLO parameters. From Fig. 2, it is apparent 

(2) (2) 

that iPi)ra(®) reaches a plateau at large A. ^ A (0) deviates from its asymptotic value by 
an amount that is of order 1/A. It is the value of V'dra(O) at t ne plateau that corresponds 
to the dimensionally- regulated matrix element [Eq. (32)]. A fit of the LO-parameter results 

(2) 

for ?/4)r,a(0) at A = 0.6, 0.7, and 0.8 to a constant plus a term that is proportional to 1/A 
yields Vdr,a(°) = °- 123 GeV 7/2 , 0.108 GeV 7/2 , and 0.095 GeV 7/2 , respectively. A similar 
fit to the NLO-parameter results at A = 0.9, 1.0, and 1.1 yields ^dr,a(°) = °- 144 GeV 7/2 , 
0.124 GeV 7 / 2 , and 0.107 GeV 7 / 2 , respectively. 

Next we consider the direct method for calculating V'drW- Substituting the LO potential- 
model parameters into Eq. (42), we obtain Vdr(°) = °- 124 GeV 7 / 2 , 0.109 GeV 7 / 2 , and 
0.096 GeV 7 / 2 at A = 0.6, 0.7, and 0.8, respectively. The results for the NLO-parameters are 
^(0) = 0.146 GeV 7 / 2 , 0.127 GeV 7 / 2 , and 0.109 GeV 7 / 2 at A = 0.9, 1.0, and 1.1, respec- 

(2) 

tively. These results for ^bk(O) are ^ n g°°d agreement with those from the extrapolation of 
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V4)ra(0) to A = oo in the hard-cutoff method. 

We consider the NLO value for tp(0) to be slightly more reliable than the LO value. 
Therefore, we use the NLO parameters to compute our central value for V4m(0)- We use 

(2) 

the LO parameters to give an indication of the uncertainty in V4)R,(0) that arises from the 
uncertainty in ip(0). As we have already discussed in Sec. II, the potential-model values of 
the string tension that derive from the NLO parameters for A = 1.0 and 1.1 span the range 
of lattice values for the string tension. The potential-model values of a s that derive from 
the NLO parameters for A = 1.0 and 1.1 are somewhat larger than the lattice values for the 
fixed parameter a s , but are compatible with the value of the running a s at the scale of the 
heavy-quark momentum m c v . We consider the lattice value of the string tension to be more 
relevant than the lattice value of the fixed parameter a s , since the latter does not take into 
account the running of a s in QCD. Therefore, we determine the central value of V'drW by 
taking the average of the values of ^1^(0) for the NLO parameters at A = 1.0 and 1.1. We 
take the difference of these values as the uncertainty that is attributable to the uncertainty 
in the potential-model parameters. We take the difference between the average of the values 

(2) (2) 

of ^drW f° r the NLO parameters at A = 1.0 and 1.1 and the average of the values of Vt,r(0) 
for the LO parameters at A = 0.7 and 0.8 to be the uncertainty that is attributable to the 
uncertainty in ip(0). We add these uncertainties in quadrature. We include an additional 
30% uncertainty to account for the fact that our potential model neglects terms of relative 
order v 2 0.3. Using the direct-method values for V'drW' we obtain 

^(0) = 0.118 ± 0.024 ± 0.035 GeV 7/2 , (52) 

where the second uncertainty arises from the v 2 error. This is the dominant source of error. 

The effects of the uncertainty in the value of ^(0) tend to cancel in the ratio ip^(0)/i/j(0). 
We compute this ratio for the LO-parameter and NLO-parameter direct-method results 
and use the same method for determining the central value and the uncertainties that we 

(2) 

described above for t/> D r(0). The result is 

V4r(0)M0) = m\v 2 = 0.50 ± 0.09 ± 0.15 GeV 2 , (53) 

where, again, the second uncertainty arises from the v 2 error. Taking m c = 1.4 GeV in 
Eq. (53), we have (v 2 ) = 0.25 ± 0.05 ± 0.08, which is in good agreement with expectations 
from the NRQCD f-scaling rules. 

19 



We can now address the numerical importance of the two-loop contribution to A^^(O) 
that is given in Eq. (36). Taking the average of the values of 7c for A = 1.0 and 1.1 in Table II 
and using the NLO value ^(0) = 0.23629 GeV 3 / 2 , we find that A^ 2) (0) = -0.048 GeV 7 / 2 . 

(2) 

This is about —40% of the value of ^^(O) in Eq. (52). 

(2) 

The v -scaling rules of NRQCD state that the hard-cutoff matrix element ip A (0), eval- 
uated at A ~ mv, should be of order m 2 v 2 ip{Q). Taking m c = 1.4 GeV, v 2 = 0.3, and 
ip(0) = 0.23629 GeV 3/2 , we obtain m 2 c v 2 ^(0) « 0.14 GeV 7/2 , which is in reasonably good 
agreement with the value of ip^\o) at A = m c v pa m for the NLO results at A = 1.0 and 1.1 
in Fig. 2. 

(2) 

In lattice determinations of ^dra(O), the lattice ultraviolet cutoff, which is of order 71 
divided by the lattice spacing, corresponds approximately to the hard cutoff A. Existing 
lattice determinations of ^dr^O) (Ref. [3]) have been carried out in the vicinity of A = m c . 
This value of A is at the boundary of the region in which asymptotic freedom allows one to 
evaluate quantities in QCD in perturbation theory. However, as can be seen from Fig. 2, 
A = m c is far from the region in which A (0) approaches its asymptotic value. Apparently, 

at A = m c , power corrections of order m c v/A are still important. Therefore, we expect that 

(2) 

the lattice determinations of V'dr(O) contain large 1/A errors. 
IX. SUMMARY 

In this paper we have presented two methods for calculating NRQCD matrix elements 
that are proportional to ^^(0), the negative of the second derivative of the wave function 
at the origin. These matrix elements enter into the relativistic corrections for the decay 
and production of heavy-quarkonium states. The matrix elements are linearly ultraviolet 
divergent, and, hence, they must be regulated. We compute the matrix elements in di- 
mensional regularization, since that is the regularization that is used most commonly in 
phenomenology. 

One method that we have presented makes use of a hard-cutoff regulator as an intermedi- 
ate step, then employs perturbation theory to compute the difference between the hard-cutoff 
matrix element and a dimensionally-regulated matrix element. This method is quite general, 
in that it requires knowledge only of the QQ wave function in the quarkonium state. In 
principle, it could be applied to wave functions that are determined directly by lattice meth- 
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ods. It involves computing cancelling quantities in the limit in which the hard cutoff is taken 
to infinity. This method has the disadvantage that the limiting procedure and cancellation 
make it difficult to achieve adequate numerical accuracy. In the case of a pure Coulomb 
potential, we have obtained analytic expressions for the quantities that enter into the hard- 
cutoff method. The result for the matrix elements agrees with the Gremm-Kapustin relation 
in that case. 

The second method that we have presented involves a computation of the matrix elements 
directly in dimensional regularization. The method is specific to potential models. It yields 
the Gremm-Kapustin relation for the matrix elements. This is not surprising, since the 
Gremm-Kapustin relation is based on the NRQCD equations of motion at leading-order 
in v, and the potential model respects those equations of motion. This method is easily 
generalized to the computation of matrix elements that are proportional to higher powers 
of V 2 acting on the wave function at the origin. We have used the results for such matrix 
elements to write a formula that resums certain contributions of higher order in v to NRQCD 
amplitudes. 

We have used both the direct method and the hard-cutoff method to evaluate the 
dimensionally-regulated quantity Vt>r.(0) f° r the ■VV' ( or the rj c ) in the Cornell potential 
model. Since the Cornell potential model contains no spin dependence, we do not distin- 
guish between the J/ip and r\ c matrix elements. If the potential itself were exact, then the 
potential model would reproduce QCD up to corrections of relative order v 2 . Existing lattice 
measurements of the static QQ potentials yield values for the string tension. In order to 
estimate errors from our choice of potential, we have carried out computations for sets of 
Cornell-potential parameters that bracket the lattice values of the string tension. 

The two methods for computing Vt)r(O) yield results that agree well numerically. Our 
final result, including estimates of uncertainties, is given in Eq. (52). The first error estimate 
arises from the uncertainty in the potential-model parameters and from the uncertainty in 
the value of the wave function at the origin that is obtained from the leptonic width of the 
J/ip. The second error estimate accounts for a 30% uncertainty that arises from the fact 
that we have neglected corrections of relative order v 2 . This is the dominant source of error. 
The effects of the uncertainty in the wave function at the origin tend to cancel in the ratio 
of matrix elements ip^(0)/ip(0) in Eq. (53). From the ratio in Eq. (53), we estimate that 
(v 2 ) = 0.25 ± 0.05 ± 0.08, which is in good agreement with the f-scaling rules of NRQCD. 
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From our analysis, it is clear that, in the hard-cutoff method, there are large corrections 
of order m c v/A, even at cutoffs A m c . This implies that existing lattice computations of 
■0dr(O) (Ref. [3]) contain large errors that arise from such power corrections. 

We also have identified a large contribution to the difference between the dimensionally- 
regulated matrix element V'drC'-O an d the hard-cutoff matrix element ip^(0) that arises 
from a two-loop contribution to a particular short-distance coefficient [Eq. (36)]. That 
contribution is suppressed only as a 2 s /v 2 relative to ip^(0) an d is about —40% of V'dr(O) in 
the case of the J/ip. Therefore, the two-loop contribution to this short-distance coefficient 
may be important numerically when one converts a lattice-regulated value of ip( 2 \0) to a 
dimensionally-regulated value. 
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APPENDIX A: TEST OF THE METHOD WITH A COULOMB WAVE FUNC- 
TION 

(2) 

In this appendix, we apply our methods for calculating V4>r(0) to the case of the wave 
function for a pure Coulomb potential. In this case, we can evaluate all of the relevant 
expressions analytically. Ultimately, we compare our results with the prediction of the 
Gremm-Kapustin relation. 

Coulomb wave functions in coordinate space and momentum space are 

tl>c{x) = ^c(0)e^ cr , (Ala) 
MP) = J <Pxe-**M*) = 8 ^c(0) (p2 J° 7 , )2 , (Alb) 
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where 7c is defined in Eq. (12). The corresponding binding energy is 

e B = -7c/™- (A2) 

It then follows from Eq. (18) that 

7b = 7c- (A3) 

We first carry out the calculation of V4)R,(0) m the hard-cutoff method. Substituting the 
Coulomb wave function into either Eq. (13) or Eq. (15), we obtain the hard-cutoff matrix 
element ip^\o): 



Vi 2) (0)=A^ c (0) 



1 - 



A 



2 



(A4) 



(A + 7c) 2 . 
(2) 

The calculation of (0) can be carried out by using Eqs. (27) and (29a). For the difference 
between the hard-cutoff and dimensionally-regulated matrix elements, we have 

(2W_.,. ,ns 2 7C ( A 4 7 B 



AV4T(0) = ^c(0) A9 2 1 — , (A5) 

B ; A 2 - 7 2 VA + 7c 7b + 7c/ 

where we have used the definition (12) to replace the prefactor a s CFiri with 27c. In the 
limit 7b — > 0, Eq. (A5) becomes 

A^ 2 b(0) = A^ c (0)-^-. (A6) 
A + 7c 

It follows from Eq. (31) that Vdra(°) 

is given by 

^1a(0) = ^ 2) (0) - A^ 2 > (0) = - 7c ^c(0) x (^-^) 2 • (A7) 

Taking the limit A — > 00, we obtain the dimensionally-regulated matrix element ■0dr(O) : 

^(0) = lim ^1 A (0) = -7c^c(0). (A8) 

We note that the deviations from the asymptotic value go as 7c/A ~ mv/A as A — > 00. 
Incidentally, had we retained the effects of the binding energy by using Aif)^ (0) instead of 

(2) 

A^n B (0), then we would have obtained 

^Ia(o) = 4 2) (o) - A^ 2) (0) = -7§lM0), (A9) 

where we have made use of the Coulomb-potential relation 7b = 7c- In this case, the A 
dependence would have vanished. However, this simplification is a special property of the 
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pure Coulomb case that arises from the fact that the Coulomb-gluon exchange in Aip^ (0) 
is, in this special case, an interaction of the complete potential. 

The result in Eq. (A8) agrees with the result from the direct method in Eq. (42) and 
with the Gremm-Kapustin relation (43). We note that, had we failed to include the two- 
loop correction to A^^(0) in Eq. (36), then we would have obtained V'dr(O) = ~~ ^Ic^ci^), 
which does not agree with the Gremm-Kapustin relation. 

APPENDIX B: INTEGRALS 

Here we record some integrals that are useful in deriving the expressions in this paper. 
1. Loop integrals 



d 3 k 1 % . /A — %\p\ 

log 



(2ti) 3 (k-p) 2 (k 2 + A 2 + ie) 8tt|p| * \A + i\p 



An\p\ \ A 

d 3 k 1 1 



arctan ( ) , (Bla) 



(2tt) 3 (k - p) 2 {k 2 + A 2 + ie) 2 8ttA(p 2 + A 2 ) ' 
d 3 k arctan (|fc| /A) 1 



(27r) 3 |fc|(fc 2 + A 2 + k) 2 16ttA 2 ' 
Analytically continuing A to iA, we obtain the following formulas: 



/ 



d 3 k 1 1 

arctan 



(2n) 3 (k - p) 2 (k 2 -A 2 + ie) 4n\p\ ' ' \iA 



d 3 k 



An\p\ \ A 

—% 



(2tt) 3 [k - p) 2 ik 2 - A 2 + ie) 2 8ttA(p 2 - A 2 ) ' 
d 3 k arctan (—i\k\/A) 1 
(2tt) 3 |fc|(fc 2 — A 2 + ie) 2 = _ 16ttA 2 ' 



(Bib) 
(Blc) 



tanh" 1 )^), (B2a) 



(B2b) 
(B2c) 
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2. Fourier Transformation 



f d 3 p arctan(|p|/A) e~ Ar 

rf3p e ^ 1 = (B3b) 



(2tt) 3 (p 2 + A 2 ) 47rr 
d 3 P i V - x 1 e 



-Ar 



^"tZTTT™ = oZT' (B3c) 



(2tt) 3 (p 2 + A 2 ) 2 8ttA' 



where r = Iccl. 
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